// header file that defines the elements and mesh struct
# pragma once

#include<array>
#include<vector>
#include<map>
#include<istream>
#include<ostream>
#include<memory>
#include<string>
#include<Eigen/Dense>
using namespace std;

struct Point {
    array<double, 3> coords; // stores (x, y, z) coordinates of a point
    int pointType; // stores the type of point (-1 means extended)
    friend istream& operator>>(istream& is, Point& p);
	friend ostream& operator<<(ostream& os, const Point& p);
};

struct Vector3D {
    double x, y, z; // components
    Vector3D() {}
    Vector3D(double x, double y, double z) : x(x), y(y), z(z) {}
    Vector3D(const Point& p1, const Point& p2) : x(p2.coords[0] - p1.coords[0]), y(p2.coords[1] - p1.coords[1]), z(p2.coords[2] - p1.coords[2]) {}

    Vector3D cross(const Vector3D& v) const;
    double inner(const Vector3D& v) const;
};

// Base struct for elements
struct Element {
public:
	map<string, double> data; // stores data value for each element, string is the name of attribute
    virtual ~Element() = default;
    int rockType; // 1 means reservoir
};

// struct element that stores nodes info
struct HexElement : public Element {
    array<size_t, 8> nodes_id; // indices for nodes -- 0, 1, 2, 3 forms a plane counterclockwise, and 4, 5, 6, 7 is above/below 0, 1, 2, 3
    array<double, 3> coords_lo; // lowest coords
    array<double, 3> coords_hi; // highest coords
    array<Vector3D, 6> normals_inward; // inward normal vectors of 6 faces
    array<Eigen::Vector3d, 6> faceCentroids; // face centroids of 6 faces

	static constexpr array<array<size_t, 4>, 6> faces = { { {0, 1, 2, 3}, {1, 5, 6, 2}, {3, 2, 6, 7}, {4, 7, 6, 5}, {0, 3, 7, 4}, {0, 4, 5, 1} } };
    // faces[i][j] -- face i node j, arranged in an order such that the right-handed normal points inward
    static constexpr array<size_t, 6> start_idx = { 0, 1, 3, 4, 0, 0 }; // starting node index for the normal vector
};

// struct element that stores nodes info for tetrahedral elements
struct TetElement : public Element {
    array<size_t, 4> nodes_id; // indices for nodes
	array<double, 3> center; // center of the element
};

struct Mesh {
    vector<Point> nodes;
    vector<unique_ptr<Element>> elems;
	virtual void initialize() = 0;
    void readVTK(const std::string& filename);
};

struct HexMesh : public Mesh {
    size_t Nx, Ny, Nz; // number of cells in each dimension
public:
    void initialize() override; // compute coords_lo and coords_hi, and compute normal vectors of all elements
    bool is_in_elem(const Point& p, const HexElement& elem) const; // check whether point is inside HexElement
    int find_containing_cell(const Point& p) const; // find the cell that contains the point
    double interpolate(const Point& p, const string& s) const; // use interpolation to find data value at point p
private:
// 使用四个节点ID作为key来存储面法向量
    map<array<size_t,4>, Eigen::Vector3d> face_normals_cache;
    static constexpr double TOLERANCE = 1e1;
};

struct TetMesh : public Mesh {
	void initialize() override; // compute element centers
};

int mesh();